home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 2002 November / SGI IRIX Base Documentation 2002 November.iso / usr / share / catman / p_man / cat3 / SCSL / clatbs.z / clatbs
Encoding:
Text File  |  2002-10-03  |  8.8 KB  |  265 lines
  1.  
  2.  
  3.  
  4. CCCCLLLLAAAATTTTBBBBSSSS((((3333SSSS))))                                                          CCCCLLLLAAAATTTTBBBBSSSS((((3333SSSS))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      CLATBS - solve one of the triangular systems  A * x = s*b, A**T * x =
  10.      s*b, or A**H * x = s*b,
  11.  
  12. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  13.      SUBROUTINE CLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE,
  14.                         CNORM, INFO )
  15.  
  16.          CHARACTER      DIAG, NORMIN, TRANS, UPLO
  17.  
  18.          INTEGER        INFO, KD, LDAB, N
  19.  
  20.          REAL           SCALE
  21.  
  22.          REAL           CNORM( * )
  23.  
  24.          COMPLEX        AB( LDAB, * ), X( * )
  25.  
  26. IIIIMMMMPPPPLLLLEEEEMMMMEEEENNNNTTTTAAAATTTTIIIIOOOONNNN
  27.      These routines are part of the SCSL Scientific Library and can be loaded
  28.      using either the -lscs or the -lscs_mp option.  The -lscs_mp option
  29.      directs the linker to use the multi-processor version of the library.
  30.  
  31.      When linking to SCSL with -lscs or -lscs_mp, the default integer size is
  32.      4 bytes (32 bits). Another version of SCSL is available in which integers
  33.      are 8 bytes (64 bits).  This version allows the user access to larger
  34.      memory sizes and helps when porting legacy Cray codes.  It can be loaded
  35.      by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
  36.      only one of the two versions; 4-byte integer and 8-byte integer library
  37.      calls cannot be mixed.
  38.  
  39. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  40.      CLATBS solves one of the triangular systems A * x = s*b, A**T * x = s*b,
  41.      or A**H * x = s*b, with scaling to prevent overflow, where A is an upper
  42.      or lower triangular band matrix.  Here A' denotes the transpose of A, x
  43.      and b are n-element vectors, and s is a scaling factor, usually less than
  44.      or equal to 1, chosen so that the components of x will be less than the
  45.      overflow threshold.  If the unscaled problem will not cause overflow, the
  46.      Level 2 BLAS routine CTBSV is called.  If the matrix A is singular
  47.      (A(j,j) = 0 for some j), then s is set to 0 and a non-trivial solution to
  48.      A*x = 0 is returned.
  49.  
  50.  
  51. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  52.      UPLO    (input) CHARACTER*1
  53.              Specifies whether the matrix A is upper or lower triangular.  =
  54.              'U':  Upper triangular
  55.              = 'L':  Lower triangular
  56.  
  57.      TRANS   (input) CHARACTER*1
  58.              Specifies the operation applied to A.  = 'N':  Solve A * x = s*b
  59.              (No transpose)
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. CCCCLLLLAAAATTTTBBBBSSSS((((3333SSSS))))                                                          CCCCLLLLAAAATTTTBBBBSSSS((((3333SSSS))))
  71.  
  72.  
  73.  
  74.              = 'T':  Solve A**T * x = s*b  (Transpose)
  75.              = 'C':  Solve A**H * x = s*b  (Conjugate transpose)
  76.  
  77.      DIAG    (input) CHARACTER*1
  78.              Specifies whether or not the matrix A is unit triangular.  = 'N':
  79.              Non-unit triangular
  80.              = 'U':  Unit triangular
  81.  
  82.      NORMIN  (input) CHARACTER*1
  83.              Specifies whether CNORM has been set or not.  = 'Y':  CNORM
  84.              contains the column norms on entry
  85.              = 'N':  CNORM is not set on entry.  On exit, the norms will be
  86.              computed and stored in CNORM.
  87.  
  88.      N       (input) INTEGER
  89.              The order of the matrix A.  N >= 0.
  90.  
  91.      KD      (input) INTEGER
  92.              The number of subdiagonals or superdiagonals in the triangular
  93.              matrix A.  KD >= 0.
  94.  
  95.      AB      (input) COMPLEX array, dimension (LDAB,N)
  96.              The upper or lower triangular band matrix A, stored in the first
  97.              KD+1 rows of the array. The j-th column of A is stored in the j-
  98.              th column of the array AB as follows:  if UPLO = 'U', AB(kd+1+i-
  99.              j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j)
  100.              = A(i,j) for j<=i<=min(n,j+kd).
  101.  
  102.      LDAB    (input) INTEGER
  103.              The leading dimension of the array AB.  LDAB >= KD+1.
  104.  
  105.      X       (input/output) COMPLEX array, dimension (N)
  106.              On entry, the right hand side b of the triangular system.  On
  107.              exit, X is overwritten by the solution vector x.
  108.  
  109.      SCALE   (output) REAL
  110.              The scaling factor s for the triangular system A * x = s*b,  A**T
  111.              * x = s*b,  or  A**H * x = s*b.  If SCALE = 0, the matrix A is
  112.              singular or badly scaled, and the vector x is an exact or
  113.              approximate solution to A*x = 0.
  114.  
  115.      CNORM   (input or output) REAL array, dimension (N)
  116.  
  117.              If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains
  118.              the norm of the off-diagonal part of the j-th column of A.  If
  119.              TRANS = 'N', CNORM(j) must be greater than or equal to the
  120.              infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be
  121.              greater than or equal to the 1-norm.
  122.  
  123.              If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns
  124.              the 1-norm of the offdiagonal part of the j-th column of A.
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. CCCCLLLLAAAATTTTBBBBSSSS((((3333SSSS))))                                                          CCCCLLLLAAAATTTTBBBBSSSS((((3333SSSS))))
  137.  
  138.  
  139.  
  140.      INFO    (output) INTEGER
  141.              = 0:  successful exit
  142.              < 0:  if INFO = -k, the k-th argument had an illegal value
  143.  
  144. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  145.      A rough bound on x is computed; if that is less than overflow, CTBSV is
  146.      called, otherwise, specific code is used which checks for possible
  147.      overflow or divide-by-zero at every operation.
  148.  
  149.      A columnwise scheme is used for solving A*x = b.  The basic algorithm if
  150.      A is lower triangular is
  151.  
  152.           x[1:n] := b[1:n]
  153.           for j = 1, ..., n
  154.                x(j) := x(j) / A(j,j)
  155.                x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
  156.           end
  157.  
  158.      Define bounds on the components of x after j iterations of the loop:
  159.         M(j) = bound on x[1:j]
  160.         G(j) = bound on x[j+1:n]
  161.      Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
  162.  
  163.      Then for iteration j+1 we have
  164.         M(j+1) <= G(j) / | A(j+1,j+1) |
  165.         G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
  166.                <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
  167.  
  168.      where CNORM(j+1) is greater than or equal to the infinity-norm of column
  169.      j+1 of A, not counting the diagonal.  Hence
  170.  
  171.         G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
  172.                      1<=i<=j
  173.      and
  174.  
  175.         |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
  176.                                       1<=i< j
  177.  
  178.      Since |x(j)| <= M(j), we use the Level 2 BLAS routine CTBSV if the
  179.      reciprocal of the largest M(j), j=1,..,n, is larger than
  180.      max(underflow, 1/overflow).
  181.  
  182.      The bound on x(j) is also used to determine when a step in the columnwise
  183.      method can be performed without fear of overflow.  If the computed bound
  184.      is greater than a large constant, x is scaled to prevent overflow, but if
  185.      the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a
  186.      non-trivial solution to A*x = 0 is found.
  187.  
  188.      Similarly, a row-wise scheme is used to solve A**T *x = b  or A**H *x =
  189.      b.  The basic algorithm for A upper triangular is
  190.  
  191.           for j = 1, ..., n
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. CCCCLLLLAAAATTTTBBBBSSSS((((3333SSSS))))                                                          CCCCLLLLAAAATTTTBBBBSSSS((((3333SSSS))))
  203.  
  204.  
  205.  
  206.                x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
  207.           end
  208.  
  209.      We simultaneously compute two bounds
  210.           G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
  211.           M(j) = bound on x(i), 1<=i<=j
  212.  
  213.      The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add
  214.      the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.  Then the
  215.      bound on x(j) is
  216.  
  217.           M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
  218.  
  219.                <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
  220.                          1<=i<=j
  221.  
  222.      and we can safely call CTBSV if 1/M(n) and 1/G(n) are both greater than
  223.      max(underflow, 1/overflow).
  224.  
  225.  
  226. SSSSEEEEEEEE AAAALLLLSSSSOOOO
  227.      INTRO_LAPACK(3S), INTRO_SCSL(3S)
  228.  
  229.      This man page is available only online.
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.